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We use molecular dynamics simulations to investigate dynamic heterogeneities and the potential 
energy landscape of the Gaussian core model (GCM). Despite the nearly Gaussian statistics of parti¬ 
cles’ displacements, the GGM exhibits giant dynamic heterogeneities close to the dynamic transition 
temperature. The divergence of the four-point susceptibility is quantitatively well described by the 
inhomogeneous version of the Mode-Goupling theory. Furthermore, the potential energy landscape 
of the GGM is characterized by a geometric transition and large energy barriers, as expected from 
the lack of activated, hopping dynamics. These observations demonstrate that all major features of 
mean-field dynamic criticality can be observed in a physically sound, three-dimensional model. 
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I. INTRODUCTION 

Supercooled liquids are characterized by collective dy¬ 
namic fluctuations, known as dynamic heterogeneities, 
which occur over longer time- and length-scales as the 
glass transition temperature Tg is approached. At the 
molecular scale, these fluctuations imply the correlated 
motion of an increasingly large number of molecules as 
relaxation slows down. To quantify dynamic hetero¬ 
geneities, a general formalism based on multi-point 
namic correlations was developed over the last years [l|. 
In particular, the four-point dynamic susceptibility Xi{t) 
allows one to evaluate the amplitude of the dynamic fluc¬ 
tuations in numerical simulations 0 and, at the cost of 
some approximations, in experiments Q. 

Despite these advances, predicting the temperature 
evolution of dynamic heterogeneities remains a big chal¬ 
lenge and none of the theories proposed so far is conclu¬ 
sive, as they describe esroerimental and numerical results 
equally well or poorly [j, Q . Amongst them, the mode¬ 
coupling theory (MCT) is known to be a microscopic and 
first principles theory of the glass transition Q. MCT 
was initially formulated as a theory of caging in liquids 
and focused on two-point correlators. A recent general¬ 
ization of MCT to inhomogeneous systems (IMCT) en¬ 
ables one to evaluate multi-point correlation functions 
and make quantitative predictions for dynamic hetero¬ 
geneities Q. Within this framework, both relaxation 
times and dynamic fluctuations diverge algebraically at 
the dynamic transition temperature Tc- These diver¬ 
gences are however “avoided” in real glass-formers: The 
dynamics at low temperature is instead governed by ther¬ 
mal activation, with a distinct super-Arrhenius tempera¬ 
ture dependence. In this regime, dynamic heterogeneities 
are expected as a manifestation of cooperatively rear- 
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ranging regions, or mosaics, as predicted by the classic 
Adam-Gibbs scenario Q. 

The random first order transition theory (RFOT), 
which was originally inspired by mean-field models of 
spin glasses, integrates these two apparently distinct sce¬ 
narios d, [i3- According to RFOT, the dynamic tran¬ 
sition predicted by MCT corresponds to the trapping of 
the system in one of the basins of its rugged free energy 
landscape. In the mean-field limit, the dynamics is com¬ 
pletely frozen-in at whereas in finite dimensions the 
transition is rounded by thermal activation and becomes 
a mere crossover. Seen from this perspective, the dy¬ 
namic transition marks a change in the topology of the 
landscape (also known as “geometric transition” [HI): 
above Tc, the system mostly resides close to saddles, 
whereas below Tc it is trapped close to local minima. The 
corresponding real space picture implies the existence of 
two distinct length scales characterizing dynamic hetero¬ 
geneities: a dynamic one, which grows algebraically on 
approaching Tc, and a static one corresponding, crudely 
speaking, to the size of the mosaics. 

Ex per imental data at ambient [l^ . [T5 | and high pres¬ 
sure [ij, as well as recent simulation results [T^ hint 
at a dynamic crossover at a temperature higher than 
Tg. Around the crossover, however, other physical mech¬ 
anisms, such as dynamic facilitation [l6l | and/or local 
structure formation [l^ , may play an important role and 
compete with the mean-field scenario. Indeed, several 
predictions of the IMCT/mean-field framework remain 
so far undetected. First, the fitted power law exponents 
describing the growth of Xi ^nd the dynamic correlation 
length do not completely agree with the ones predicted 
by IMCT [l^. Moreover, according to MCT and IMCT, 
dynamic fluctuations grow near Tc but the single-particle 
dynamics remains essentially Gaussian. This counter¬ 
intuitive behavior is absent in standard glass-formers, for 
which Xi the non-Gaussian parameter 02 are typi¬ 
cally correlated and grow concomitantly as the dynam¬ 
ics slows down M- Finally, the saddles that become 
marginally stable at Tc should be delocalized M , but no 
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trace of such extended modes was detected in the poten¬ 
tial energy surface of common glass-formers [^ . 

In this paper, we put the mean-field scenario to a cru¬ 
cial test by studying the approach to the dynamic tran¬ 
sition in the Gaussian core model (GCM) [2l|- Con¬ 
ventional model glass formers, such as Lennard-Jones 
and hard sphere fluids, are characterized by short-range, 
harshly repulsive potentials. In the GCM, instead, par¬ 
ticles interact via an ultra-soft repulsive potential v{r) = 
whose tail plays a key role at high 
density. At sufficiently high density, the GCM becomes 
a good glass-former and its average dynamics is well de¬ 
scribed by MCT m, [U . Simulations of the high-density 
GCM are computationally expensive, due to the large 
number of neighbors each particle is interacting with. 
Therefore, characterizing dynamic fluctuations and the 
energy landscape of this model requires a major overhaul 
of the numerical protocol compared to previous stud¬ 
ies [H n. By employing state-of-art molecular dy¬ 
namics simulations on graphics processor units (GPU), 
we demonstrate that the GCM provides a striking in¬ 
carnation of mean-field dynamic criticality in a three- 
dimensional system and provides a solid reference to un¬ 
derstand how and when the mean-field scenario is washed 
out in more common glass-formers. 


II. METHODS 

We use molecular dynamics simulations in the NVT 
ensemble with a Nose-Hoover thermostat to study N = 
4000 monodisperse GCM particles. The potential is cut 
and shifted and smoothed at Tc = d.Str with the XPLOR 
cutoff [ 2 ^, which ensures continuity of the force at the 
cutoff. For the energy landscape analysis we choose a 
smaller system {N = 2000). In the following, we will 
use cr, 10“®e, and 10“®e/fc_B as units of the length, en¬ 
ergy, and temperature, respectively. We focus on super¬ 
cooled fluids along the isochore p = 2.0, for which the 
thermodynamically stable state is the BCG crystal (for 
T < 8.2) [2^ [2^. The crystallization kinetics is very 
slow and virtually negligible in our simulations [2^, l24j| . 
We note that the density is much higher than p k, 0.24, 
the re-entrant melting density of GCM [26l - l^ . In the 
low density limit, the GCM approaches asymptotyically 
the three-dimensional hard sphere model [^. In this 
regime, the physics of model differs markedly from the 
one observed at high density. To ensure good statistics 
on the four-point dynamic susceptibility, we performed 
four independent production runs for each temperature 
and the simulation time for each trajectory was typically 
100 times longer than the structural relaxation time Tq, 
(see below). 

Computer simulations of the GCM in the high den¬ 
sity, supercooled regime are computationally demanding 
and require particular care. Due to quasi-long range na¬ 
ture of the potential and the large statistics needed to 
evaluate X 4 , the simulation protocol must be efficient. 


Moreover, in contrast to systems with harshly repulsive 
interactions, the force summation is not dominated by 
the first coordination shell and involves a large number 
of atoms. This, in turn, rises obvious issues of numer¬ 
ical accuracy. To tackle these issues, we employed the 
HOOMD simulation package [ 2 ^ [ 2 ^, a state-of-art 
molecular dynamics code running on graphic processor 
units (GPU) with double precision arithmetics. HOOMD 
is currently one of the few simulation codes running en¬ 
tirely on GPU that allows double precision evaluation of 
both the forces and the integration step. We checked the 
results obtained with HOOMD against those of in-house 
simulation codes running on GPU over the available tem¬ 
perature range and found them to be nicely consistent. 
An initial batch of simulations was performed using the 
RUMD simulation package [3l|. RUMD implements the 
force calculation in single precision, which turns out to 
be insufficient for the GCM at high density. In fact, 
small but systematic differences between the simulations 
performed with RUMD and with our in-house codes ap¬ 
peared at sufficiently low temperature, in both NVE and 
NVT ensembles. These data sets were not retained in 
the analysis. 

The stationary points of the potential energy surface 
(PES) were located using standard numerical strategies 
adopted in earlier studies on LJ mixtures m- To locate 
local minima and saddle points we minimized the total 
potential energy U and the total force squared 

w = (1) 

i=l 

where fi is the norm of the force vector on particle i re¬ 
spectively, using the LBFGS minimization algorithm . 
For each studied temperature, we considered 80 indepen¬ 
dent configurations as starting points of our U- and W- 
minimizations. 

Due to the large system size and the long range cutoff, 
U- and W-minimizations for the GCM are technically 
more difficult than for LJ mixtures. To reduce the com¬ 
putational burden without biasing the results, we used a 
smaller system size {N = 2000) compared to the one used 
to characterize the dynamics (iV = 4000). We note that 
for N = 2000 particle at p = 2.0 the box length L = 10 
is only slightly larger than twice the cutoff distance. 

Due to softness of the potential and the large num¬ 
ber of interacting particles, minimizations require high 
numerical precision as well as a smooth cut-off scheme. 
Thanks to the XPLOR cutoff, we could converge (7- 
minimizations to values of the mean squared total force 
of order W ~ 10“^^. For most of the configurations 
we located this way, the dynamical matrix contained no 
imaginary modes, ^-minimizations that did not converge 
to true minima were therefore discarded from the analy¬ 
sis. The fraction of discarded configurations ranged from 
less than 10% (close to Tc) to about 20% (at the highest 
temperatures). We note however that inclusion of such 
spurious configurations does not alter the average inher¬ 
ent structure energy within statistical noise. To try to 
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FIG. 1: Relaxation times Ta (red filled circles) and inverse 
diffusion constants D~^ (red open circles) of the GCM as 
a function of the reduced temperature e = T/Tc — 1. The 
D~^ of the GCM obtained in Ref. [2^ (black crosses) are also 
plotted. The dashed lines are fits to Ta ~ e ^ with 7 = 2.7 
and D~^ « with 7 = 2 . 6 . 

improve the convergence, we replaced the XPLOR cut¬ 
off by a smoother, cubic interpolation scheme m dur¬ 
ing the minimization. We found that the cubic splined 
cutoff did not improve appreciably the accuracy of the 
minimizations compared to the XPLOR cutoff. This in¬ 
dicates that the main numerical difficulty lies in the high 
dimensionality of the system, which makes the minimiza¬ 
tion problem ill-conditioned. 

VP-minimizations are known to locate true saddle 
points only rarely [1^. Most of the points located 
through such a procedure are, in fact, gwasl-saddles, i.e. 
local minima of W with finite W values. Our minimiza¬ 
tion algorithm locates configurations with W of the order 
10“^^, which is close to but still larger than the thresh¬ 
old we used for local minima {W = 10“^^). We conclude 
that the points located by our VP-minimizations should 
be considered as quasi-saddles. Previous studies showed 
however that the statistical properties of quasi-saddles 
and true saddles are practically indistinguishable [^ . 
On this basis, we assumed the equivalence of these two 
types of points in our analysis. 

III. RESULTS 

Before investigating dynamic heterogeneities, we study 
the relaxation dynamics of the model close to Tc- We 
measured the time-dependent overlap function defined 
by {F(t)) = -«))i where ARi{t) is 

the displacement of the i-th particle in the time interval 
t, 0(x) is the Heavyside’s step function, and we choose 
a = 0.3, which maximizes the four-point susceptibility 
defined later. {F{t)) gives the average fraction of parti¬ 
cles which moved more than a after a time t. The relax¬ 
ation time Ta is defined by {F(t = Ta)) = e~^. Our equi¬ 


librium simulations extend down to T = 2.8, which cor¬ 
responds to a relaxation time 4 times longer than those 
accessible to previous simulations [H, . The increase 

of Ta becomes non-Arrhenius around the onset tempera¬ 
ture To = 5.0 and displays a power law behavior ^ 
at low temperature. 

Figure [T] compares the relaxation times Tq and the 
diffusion constants D obtained in this work to those of 
Ref. [ 2 ^. The data are plotted against the reduced tem¬ 
perature e = T/Tc — 1. The mode-coupling temperature 
To was determined following the standard fitting proce¬ 
dure: we fixed the exponent 7 = 2.7, obtained by solving 
the MCT equation for the GCM, and then determined 
To by linear regression of Ta against T in the range 
T < 3.2. The so obtained value of To = 2.68 is only 
20% smaller than the theoretical prediction (Tc = 3.2), 
which should be contrasted with more than 100% for the 
Kob-Andersen Lennard-Jones (KA) mixture [13, We 
point out that the MCT power-law fit works for a wider 
range of e, viz. down to e = 0.037, than in the KA mix¬ 
ture, for which deviations are already apparent around 
e = 0.1 [33 l35j|. Furthermore, the new simulations allow 
us to detect a slight discrepancy between the exponents 
that fit the relaxation times (7 = 2.7) and the inverse 
diffusion constants (7 = 2.6). In the latter case, we fixed 
the critical temperature to the one obtained from the 
relaxation times analysis and only adjusted the power law 
amplitude and exponent. These minor deviations from 
MCT predictions should be contrasted to those found in 
other models, which are typically around 20-25%. We 
note that the small discrepancy between r and 1/D ob¬ 
served in the GCM could be explained even within the 
framework of the MCT itself, as discussed in Ref. [s^. 

To characterize dynamic heterogeneities in the 
GCM, we consider two different observables. First, 
we evaluate the non-Gaussian parameter a 2 {t) = 

3(Ai?(t)^)/5(Ai?(t)^)^ — 1, which quantifies how much 
the particles’ displacements deviates from a Gaussian 
distribution. In the left panel of Fig. [51 we plot 02 (t) 
for the KA mixture and the GCM for similar relaxation 
time windows. a 2 {t) of the KA mixture grows rapidly 
as the temperature is decreased whereas the growth of 
a 2 {t) of the GCM is moderate. This observation is 
consistent with the shape of the van Hove functions 
Gs{r,t) = '^)) near t = Ta- Fig¬ 

ure [3] shows the probability distribution of the logarithm 
of single-particle displacements, which is related to the 
van Hove functions as P(log^Qr, t) = (In 10)47rr^Gs(r, t), 
of the GCM at T = 3.0 and 2.8. For standard model glass 
formers, the shape of this distribution function is bimodal 
close to Tj,, which proves the existence of hopping-like 
motion [^. For the GCM, by contrast, the distribution 
functions remains unimodal even at the lowest investi¬ 
gated temperature. The absence of hopping-like motion 
in the GCM is consistent with our analysis of the poten¬ 
tial energy landscape (see below) which indicates that 
thermally activated relaxation is strongly suppressed. 
We conclude that the distribution of single-particle mo- 





4 



10-'' 10° 10^ 10^ 10° lO'* 10-'' 10“ 10'' 10° 10° 10“* 


t It 

o 


t It 

o 


FIG. 2: Left: non-Gaussian parameter 02 ( 1 ) of the GGM (red 
circles) at T = 5.0, 3.4, 3.0 and 2.9, and of type A particles 
of the KA mixture (empty squares) at T = 1.0, 0.6, 0.5 
and 0.466. Right: four-point susceptibility X4(f) of fh® GGM 
(red circles) at T = 4.0, 3.0, 2.9 and 2.8, and of the KA fl^ 
(empty squares) at T = 0.55, 0.5, 0.47 and 0.45. Downward 
arrows indicate the corresponding relaxation times Ta of the 
GGM. Time t is scaled by To, the relaxation time at the onset 
temperature To. [To = 5.0 and To = 3100 for the GGM, and 
To = 1.0 and To = 15 for the KA in the unit of Ref. [33| 1 
All data sets in the right panel are obtained from molecular 
dynamics simulations in the NVT ensemble. 


bility in the GCM is homogeneous, even when the slow 
dynamics is well developed. The shape of Gg {r, t) of the 
KA mixture is bimodal, indicating coexistence of highly 
mobile and immobile particles |35l |. while that of the 
GCM remains unimodal even at the lowest temperature. 

Next, we evaluate the four-point dynamic susceptibil¬ 
ity defined by X 4 (t) = N[{F(t)'^) — (F(t))^], which nuan- 
tifies the cooperative motion of particles in fluids 
Strikingly, the trend is now reversed as shown in the right 
panel of Fig.[2J X 4 (t) grows far more strongly in the GCM 
than in the KA mixture. Around T^, the maximum of 
Xiit), called herein X 4 , in the GCM is one order of magni¬ 
tude larger than that of the KA mixture. Note that while 
the dynamic susceptibility depends in general on the sta¬ 
tistical ensemble, this result holds also in the ensemble 
where all global variables are allowed to fluctuate (see 
Appendix). Thus dynamic fluctuations in the GCM are 
significantly larger than those in other standard models 
around T^. 

The opposite trends of 02 (f) and X 4 (f) may look con¬ 
tradicting at a first glance. This implies that the nature 
of dynamic heterogeneities of the GCM is qualitatively 
different from other conventional glass formers. Large 
values of 0:2 (i) are a direct consequence of the large dis¬ 
placements of individual mobile particles and do not nec¬ 
essarily reflect the extent of cooperative motion. On the 
other hand, X 4(0 is defined as the variance of the over¬ 
lap function, which does not depend on how far the mo¬ 
bile particles have moved, and it is thus more sensitive 



FIG. 3: Probability distribution of the logarithm of the par¬ 
ticle displacements P(logjQr, t) = (In 10)47rr°Gs(r, t). Top: 
results for T = 3.0. From left to right, I/to = 4, 17, 68, 271 
and 1082, where ToijTo = 75. Bottom: results for T = 2.8. 
From left to right, I/to = 34, 135, 406, 812, 3112 and 8659, 
where Tc/to = 813. 


to the net cooperativity. Therefore, the suppression of 
a 2 {t) and the concomitant enhancement of Xi{^) m the 
GCM imply slight but spatially extended modulations of 
the mobility field. 

These inferences are confirmed by visual inspection of 
the mobility field close to the dynamical critical temper¬ 
ature. To visulalize the giant dynamic fluctuations that 
give rise to the increase of X 4 close to the dynamic transi¬ 
tion, we evaluate the mobility of the particles after time 
t as Sri{t,to) = \ri{t + to) — ri{to)\. In Figure |4] we show 
typical snapshots of the mobility field at T = 2.9, close to 
the dynamic transition. At this temperature, the max¬ 
imum of the dynamic susceptibility has reached about 
200. The radii of the spheres are proportional to the 
particles’ displacements after a time t = 4 x 10® si Ta. 
We clearly see that the mobility field is characterized by 
extended regions of either mobile or immobile particles, 
with very smooth variations over space. Note that no 
coarse-graining or time averaging is involved in our cal¬ 
culation of Sr. Visual inspection suggests that the corre¬ 
lation length is comparable to the system size at this tem¬ 
perature. This is in turn compatible with the existence 
of finite size effects around and below this temperature. 

The coexistence of giant dynamic heterogeneities and 
Gaussian-like dynamics is a strong evidence of mean- 
field dynamic criticality. According to MCT and IMGT, 
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FIG. 4: Typical snapshots of the particles’ mobility held at a temperature T = 2.9. Particles are shown as spheres of radius 
proportional to the mobility Sr{t,to) after a time t — A x 10® « Tq. The time origins to at which the conhgurations are shown 
are separated by at least 10 structural relaxation times. 


the amplitude of a 2 {t) remain small [s^, whereas xt 
diverges on approaching with an exponent that de¬ 
pends on both microscopic dynamics and statistical en¬ 
semble iam. For molecular dynamics simulations 
in the NVT ensemble, the susceptibility is predicted to 
diverge as xl Figure [5] shows the dependence of 

xt on e. The data for the KA mixture are not well de¬ 
scribed by IMCT, whereas for the GCM xl follows the 
power law predicted by the theory over the range of tem¬ 
perature 2.8 < T < 3.4. We emphasize that in Fig. [5] 
the critical temperatures Tc are those obtained by fitting 
the relaxation times data. Our data demonstrate that 
T, 1/D, and xl can be fitted by power laws over com¬ 
parable £ ranges and using the appropriate set of MCT 
exponents. 

Deviations from the IMCT prediction are observed 
only at the lowest temperature in Fig.jSJ This deviation is 
most likely due to finite size effects, which will naturally 


appear if Xii't) has a genuine divergence, rather than to 
the crossover between MCT and activated regimes. In 
the inset of Fig. [5l we plot the histogram of the overlap 
F{t) for the CCM at t = Tq,. By definition, the mean 
value of the histogram is {F{t = Tq,)) = e~^. Indeed, the 
histogram at T = 3.0 is unimodal with the mean value of 
e“^, whereas it becomes very broad at T = 2.8, suggest¬ 
ing that a correlation length becomes comparable to the 
system size. We note that such stron g fin ite size effects 
were not observed in the KA mixture |18l . l40j| . 

We now proceed to discuss the dynamics of the CCM 
from the perspective of the energy landscape [4l|,|43. Ac¬ 
cording to the RFOT scenario, the MCT crossover should 
be accompanied by a “geometric” transition in the en¬ 
ergy landscape Q: in mean-field, the unstable modes 
separating the free energy minima become marginally 
stable as Tc is approached and eventually disappear be¬ 
low a certain threshold energy m- In finite dimensions. 
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FIG. 5: Maximum value of the dynamic susceptibility, xh 
against the reduced temperature e = T/Tc — 1 in the GCM 
(red circles) and in the KA mixture |1^ (empty squares), 
where Tc = 2.68 for the GCM and Tc = 0.435 for the KA. The 
solid line is the IMCT prediction, xX ~ ■ The inset shows 

the histogram of the overlap between two conhgurations F{t) 
with the time interval t = Tc for the GCM at T = 3.0 and 
2 . 8 . 

remnants of this geometrical transition should be visi¬ 
ble in the potential energy landscape around the MCT 
crossover: Above Tc the system relaxes mostly via un¬ 
stable soft modes, while activated relaxation over energy 
barriers takes over below Tc- The crossover between the 
two regimes takes place at an energy threshold Cth, below 
which the system resides mostly close to local minima of 
the potential energy instead of stationary points of arbi¬ 
trary order. Early numerical simulations of the KA mix¬ 
ture gave support to this scenario and found a crossover 
temperature Tth very close to Tc [HI, , but were later 
on called into question [H) IHl • 

To determine the statistics of stationary points, we 
located minima and saddles by applying the LBFGS 
minimization algorithm to the total potential energy U 
and to the total force W respectively (see Sec. HU). 
Although local minima of W do not necessarily corre¬ 
spond to true stationary points [43| . IT-minimizations 
yield a fairly robust measurement of the energy thresh¬ 
old Cth and of the typical energy barriers [i^. FigurelSKa) 
shows the saddle point energy, Cgad, as obtained from W- 
minimizations, against the fraction of unstable modes / 
found in the spectrum of the dynamical matrix [^ . As 
in other models, these two quantifies are roughly linearly 
related. From a linear fit Csad = e-th + (3Ae)/ we extract 
the threshold energy Cth = 1.70. Comparing this with 
the temperature dependence of the energy of minima Cm 
(Fig.|n](b)), we obtain a threshold temperature Tth ~ 2.7 
in excellent agreement with Tc- The slope Ae gives an 
estimate of the average barrier height, since it is the en¬ 
ergy cost to increase by one the order of a stationary 
point [ll|. We found AejTc = 15.2 in the GCM, which is 
appreciably larger AejTc ~ 10 observed in other model 
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FIG. 6: Analysis of the potential energy landscape of the 
GGM. (a) Energy of saddle points Cg^d a function of the 
fraction of unstable modes /. The solid line is a linear fit 
Csad = eth + (3Ae)/ with eth = 1-70 and AejTc = 15.2. 
(b) The inherent structure energy ds vs. temperature. The 
horizontal dashed line indicates eth = 1.70. Vertical arrows 
indicate the MGT temperature Tc = 2.68 and the onset tem¬ 
perature To = 5.0. (c) Distributions of the participation ratio 
of unstable modes of the KA mixture {N — 2000) at T = 0.45 
and of the GCM at T = 2.9. 

glass formers i Ei- We estimate that the increased 
barrier height hampers the activated relaxation by a fac¬ 
tor exp(15)/exp(lO) ~ 150. As the activated relaxation 
channels are strongly suppressed, the MCT-like critical 
behavior dominates the slow dynamics of the model [HI- 
Further support for the geometric transition scenario is 
provided by the analysis of the mode localization. It has 
been argued that the unstable directions that disappear 
at the dynamic transition are delocalized in the mean- 
field scenario [l^. This contrasts with the typical ob¬ 
servation that unstable modes of common glass-formers 
become increasingly localized as T decreases [43|. To 
evaluate the spatial localization of the unstable modes, 
we calculate their participation ratio p on a per-mode 
basis 

p(w) = ^ ^(ei(w) • ei(a;))2^ ^ (2) 

where ei{uj) is the displacement of particle i along the 
eigenvector corresponding to the eigenfrequency w. In 
Fig. [HKc) we plot the distribution of the participation ra¬ 
tio of unstable modes around Tc- As in various conven¬ 
tional glass formers [l^, the unstable modes of the KA 
mixture have small participation ratios and are there¬ 
fore spatially localized. By contrast, the distribution is 
broader in the GCM and considerably shifted towards 
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larger value of p (0.6 for plane waves, 1 for completely 
delocalized modes). Visual inspection suggests that these 
modes have a complex spatial structure, in which some 
groups of particles undergo cooperative motions, while 
others display incoherent displacements. While deter¬ 
mining the precise nature of the modes in the GCM would 

While determining the precise nature of the modes in 
the GCM would require an analysis of iV-dependence of 
the spectrum, and thus of the mobility edge [i^, our 
analysis suggests that the giant dynamic heterogeneities 
of the GCM might indeed be associated to these extended 
unstable modes. In the KA mixture, instead, dynamic 
heterogeneities build up through dynamic facilitation of 
localized elementary rearrangements 0, which might 
be related to modes localized outside locally stable do¬ 
mains [13 . 

IV. CONCLUSIONS 

In summary, we presented evidence that the mean-field 
scenario predicting diverging dynamic fluctuations and a 
geometric transition close to the dynamic transition can 
be observed in a three-dimensional model system at suf¬ 
ficiently high density. Our results shed new light on the 
physical nature of dynamic heterogeneities predicted by 
IMCT. The approach to the dynamic transition is accom¬ 
panied by rapidly growing dynamic correlations and by 
the disappearance of extended unstable modes of the po¬ 
tential energy landscape. We identify a clear fingerprint 
of mean-field dynamic criticality, that is, the coexistence 
of the strong dynamic fluctuations and nearly Gaussian 
distribution of a single particle displacement. Which 
finite-dimensional features (e.g. locally preferred struc¬ 
tures, dynamic facilitation) mask the mean-field physics 
in actual supercooled liquids is a question that need to 
be addressed in future studies. 


Acknowledgments 

We thank two anonymous reviewers for useful re¬ 
marks and constructive criticisms on a previous ver¬ 
sion of the manuscript. We acknowledge the HPGLR 
Genter of Gompetence in High-Performance Comput¬ 
ing of Languedoc-Roussillon (France) and ACCMS of 
Kyoto University (Japan) for allocation of CPU time. 
A.I. acknowledges the financial support from JSPS post¬ 
doctoral fellowship for research abroad and KAKENHI 
No. 26887021. K.M. is supported by KAKENHI No. 
24340098, 25103005, and 25000002. 

Appendix A: Checks on ensemble dependence 

Most of our simulations were performed in the NVT 
ensemble by means of the Nose-Hoover thermostat. It is 
important therefore to make sure that the choice of the 
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FIG. 7: Effect of the thermostat relaxation time Ttherm on 
the dynamic susceptibility y 4 at T = 3.4 in the GCM. The 
values of the thermostat relaxation time are indicated in the 
figure caption. Most of our simulations were carried out using 
rtherm ~ 200. 
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FIG. 8: Statistical ensemble dependence of the dynamic sus¬ 
ceptibility in the GCM. The dynamic susceptibility measured 
in the NVT ensemble , circles) and in the NVE en¬ 
semble 1 filled squares) are plotted for T = 3.4. Open 

squares indicate Xi^^ + T^Xt/cv, which must be identical 
to 


thermostat relaxation time did not affect the dynamic 
quantities, in particular the dynamic susceptibility. In 
Fig. [7] we show evaluated using four different ther¬ 

mostat relaxation times Ttherm = 2, 20, 200, 2000. In our 
simulations we choose Ttherm = 200. The differences be¬ 
tween Xii^) obtained with different Ttherm remain within 
statistical uncertainties, thus our choice of Ttherm does 
not affect the results. 

As is well established, the dynamic susceptibility de¬ 
pends in general on the statistical ensemble Q. Since 
the thermostat relaxation time Ttherm is long enough such 











FIG. 9: Log-log plot of x-ii^) in th® NVT ensemble for the 
GCM (red circles) at T — 4.0, 3.0, 2.9 and 2.8 (from left to 
right) and for the KA mixture (empty squares) at T = 0.55, 
0.5, 0.47 and 0.45 (from left to right). Blue diamonds indicate 
the full dynamic susceptibility X 4 °*^* of the KA model at T = 
0.55, 0.5 and 0.47 (from left to right). 

that X 4 (t) is independent from rtherm itself (Fig. [7]), the 
evaluated X 4 (t) can be safely considered as the dynamic 
susceptibility in the NVT ensemble [i^. It is possible 
to transform the dynamic susceptibility between ensem¬ 
bles by means of well-known transformation formulas. 
In particular, the dynamic susceptibility in the NVT en¬ 
semble is related to the one in the NVE ensemble by 
+ T'^Xt/^v, where xt = d{F{t))/dT and 
cy is the specific heat at constant volume. This, in turn, 
provides us with an internal check of our calculations. In 
Fig.|8]we compare xf^^y Xa'^'^ ^ and Xi^^ FT'^XtIcv- 
For this calculation, we averaged over 16 independent 
simulation runs for each ensemble to improve the statis¬ 
tical accuracy. We see that agrees with 

within error bars. We also note that the difference 
between Xi^^ and Xi^'^ is small at this temperature. 

Finally, we point out that Xa'^'^ obtained from our 
NVT simulation is essentially the same as the total dy¬ 
namic susceptibility X 4 °*^'^' j defined as the volume integral 
of the four point correlator. X 4 °*^'^' can be expressed as 
the sum of Xi^’^ and X4|i5Afi t^e fluctuation of dynamics 
originating from particle number fluctuations. 

To assess the role of xaIsn quantitatively, we calculated 
it explicitly for both GCM and KA models. For a binary 
mixture, Xi\sN can be expressed as 

Xi\sN = X^Hi + XpXcH 2 + Xc^3 + {F{t)fH4 

+ {F{t))xpH5 + {F{t))xcHe (Al) 

with 

Hi = p^cSii + 2pVc(l - c)Sl 2 + p\l - c)S 22 

H 2 = 2pc{l — c)5'ii -I- 2p(l — 2c)\/c(l — c)S'i 2 
-2pc(I - c)S 22 

H 3 = c(I - c)^S'ii - 2c(I - c)\/c(l - c)S'i 2 
-bc^(l - c)S 22 

Hi = cSii + 2y^c(l - c)5'i2 + (1 - 0)822 

H 3 = 2^0511 -I- 4p-^c(l - c)5'i 2 + 2p(I - 0)822 

He = 2c(I - c) 8 ii + 2(1 - 2 c)i/c{l - 0)812 
—2c(l — c)8o'} 


where Xp = i9{F{t))/dp)c, Xc = {.d{F{t))/dc)p, p is the 
dimensionless number density and c is the concentration 
of the species 1 in the binary mixture. 8 nm is the A: —>■ 0 
limit of the partial structure factor 8 nm{k). The present 
expressions differ from those in Ref. , because we treat 
p (instead of the packing fraction) and c as independent 
variables. We checked that both formulations give the 
same results for X4|5Af in the model at T = 0.55. 
Finally, if one sets c = 1 in the above expressions, one 
gets the expression for X 4 |i 5 Af for monodisperse systems. 


First we calculated X4|(5Af for the GCM at T = 3.2 and 
3.4 using the expressions above and found it to be of or¬ 
der 10“^ at these temperatures and thus negligibly small. 
This is largely due to the extremely small compressibility 
(to which xa\sn is proportional, as shown in Eq. dUD). 
Indeed, we found S'(A: —>• 0) « 10 ® at the temperatures 
we studied [ 2 ^. Furthermore, Xp in Eq- m does not 
diverge faster than Xt- This can be rationalized by the 
fact that Tc scales like the melting temperature as a func¬ 
tion of density, i.e. Tc exp(—6p^/^) and thus Xp should 
grow with the same exponent as xt- From these facts, 
we conclude that X 4 °*^* ~ Xnvt for the GCM. 


Next, let us assess xaIsn for the KA model. We per¬ 
formed additional molecular dynamics simulations of the 
KA model at T = 0.55, 0.5 and 0.47 to evaluate X4|<5Af- 
The system size of these simulations was N = 1000. We 
also performed a simulation with a larger system size 
{N = 8000) at r = 0.55 to confirm our extrapolation 
of 8 nm{k) in the A: ^ 0 limit. The maximum values of 
X 4 | 5 Ar are 2.5, 7.3, and 15.9 at T = 0.55, 0.5 and 0.47, 
respectively. In Fig. [9] we plot the full dynamic suscepti¬ 
bility X 4 °*^* = Xi^'^ + X4|<5Af for the KA model for these 
three temperatures and compare them to the results in 
Fig. Oof the main text. We point out again that 
should be indistinguishable from Xi'^'^ for the GCM. 


From these observations, we conclude that the contrast 
between the two models around Tc remains sharp even 
by including the Xa\sn term and confirms that dynamic 
fluctuations are much more pronounced in the GCM than 
in the KA model. 
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